### RPS Public Opinion Experiment Analysis 
### May 16 2017

# Set up
## Set Working Directory to folder where you downloaded replication data

library(ggplot2)
library(foreign) 
library(scales) # for ggplot2 %
library(RColorBrewer) # colors
library(survey) # survey package
require(MASS)
require(Hmisc)
require(reshape2)
library(stargazer)

rm(list=ls()) #to remove all master from workspace

options(digits=10)
options(scipen=10)
set.seed(02139)

## Load cleaned survey data
load(file="rps_experiment.RData")
states<-read.csv(file="StateCodes.csv")
master<-merge(master, states, by.x="state", by.y="Name", all.x=T, all.y=T)
#####
### Main Analysis, using 4-point scale #####
#####
# OLS model without covariates
reg <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+ factor(Randomization4b) + factor(Randomization5b), data=master) #+factor(ideo7)
summary(reg)


# Ordered probit model without covariates
reg2 <- polr(factor(billsup)~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+ factor(Randomization4b) + factor(Randomization5b), data=master, Hess=TRUE, method="probit")

# OLS model with covariates, Supplementary Table 3
reg.covariates <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+ factor(Randomization4b) + factor(Randomization5b)+factor(pid3_leaners)+factor(educ)+factor(gender)+age+factor(race)+income2, data=master) 
summary(reg.covariates)

### Create weights based on census targets
age.dist <- data.frame(age2 = c("young", "middle", "old"),
                       Freq = nrow(master) * c(36.5/76,26.4/76, 13/76))

sex.dist <- data.frame(gender = c("male", "female"),
                       Freq = nrow(master) * c(0.49, 0.51))
race.dist <- data.frame(race = c("white", "black","hispanic", "other"),
                       Freq = nrow(master) * c(0.62, 0.13,0.18, 0.07))

edu.dist <- data.frame(educ = c( "Did not graduate from high school","High school graduate",  "Some college, but no degree","2-year college degree", "4-year college degree","Postgraduate degree (MA, MBA, MD, JD, PhD, etc.)"),
                       Freq = nrow(master) * c(0.117039137,	0.289542374,	0.190966762,	0.095585502,	0.194920897,	0.111937159))
master2<-subset(master, !is.na(gender) & !is.na(educ) & !is.na(race)& !is.na(age2))
master.unweighted <- svydesign(ids=~1, data=master2)
master.rake <- rake(design = master.unweighted,
                   sample.margins = list(~gender, ~educ, ~race, ~age2),
                   population.margins = list(sex.dist, edu.dist, race.dist, age.dist))
summary(weights(master.rake))
library(dplyr)
master.wtd <- master2 %>%
    mutate(weight = weights(master.rake))
    master.wtd$weight[ master.wtd$weight>4]<-4
## raw data
prop.table(svytable(~gender, design=svydesign(~1, weights=~1, data=master.wtd)))
prop.table(svytable(~educ, design=svydesign(~1, weights=~1, data=master.wtd)))
prop.table(svytable(~race, design=svydesign(~1, weights=~1, data=master.wtd)))
prop.table(svytable(~age2, design=svydesign(~1, weights=~1, data=master.wtd)))

## weighted data
prop.table(svytable(~gender, design=svydesign(~1, weights=~weight, data=master.wtd)))
prop.table(svytable(~educ, design=svydesign(~1, weights=~weight, data=master.wtd)))
prop.table(svytable(~race, design=svydesign(~1, weights=~weight, data=master.wtd)))
prop.table(svytable(~age2, design=svydesign(~1, weights=~weight, data=master.wtd)))
summary(master.wtd$weight)

# OLS model with weights
reg.weights <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+ factor(Randomization4b) + factor(Randomization5b), data=master.wtd, weight=weight) 
summary(reg.weights)

##  Supplementary Table 3
## comparing models without covariates/weights, with covariates, with weights
stargazer(reg,reg.covariates, reg.weights)

##  Supplementary Table 4
## comparing models from ols and ordered probit
stargazer(reg,reg2)

### Figures 2 in main paper
## Plot
variables<-read.csv(file="variables.csv")
variables$coef<-NA
variables$SE<-NA

#messenger
variables$coef[24]<-0
variables$coef[25]<-coef(reg)[2]
variables$coef[26]<-coef(reg)[3]

#cost
variables$coef[3]<-0
variables$coef[4]<-coef(reg)[4]
variables$coef[5]<-coef(reg)[5]
g
#jobs
variables$coef[8]<-0
variables$coef[9]<-coef(reg)[6]
variables$coef[10]<-coef(reg)[7]

#air
variables$coef[13]<-0
variables$coef[14]<-coef(reg)[8]

#climate
variables$coef[18]<-0
variables$coef[19]<-coef(reg)[9]
variables$coef[20]<-coef(reg)[10]
variables$coef[21]<-coef(reg)[11]

SE<-summary(reg)$coefficients[, 2]
#messenger
variables$SE[24]<-0
variables$SE[25]<-SE[2]
variables$SE[26]<-SE[3]

#cost
variables$SE[3]<-0
variables$SE[4]<-SE[4]
variables$SE[5]<-SE[5]

#jobs
variables$SE[8]<-0
variables$SE[9]<-SE[6]
variables$SE[10]<-SE[7]

#air
variables$SE[13]<-0
variables$SE[14]<-SE[8]

#climate
variables$SE[18]<-0
variables$SE[19]<-SE[9]
variables$SE[20]<-SE[10]
variables$SE[21]<-SE[11]

variables1<-variables
variables1$scale<-"Bill Support on Four-Point Scale"
fontface<-rev(variables$fontface[1:26])

order<-1:26
variables <- transform(variables, Categories=reorder(Categories, -order) ) 

variables2<-variables[2:14,]

title<-paste("Effects of Bill Attributes on RPS Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_design",".pdf",sep=""), width=6.33, height=5.25) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("4-point Scale for Bill Support")+
  theme(axis.title.y = element_blank())+
  theme(axis.text.y=element_text(size=10, face=rev(variables2$fontface)))  +
    scale_x_continuous(limits = c(-.5, .5),breaks = c(-.4,-.2,0,.2,.4))+
 coord_cartesian(xlim=c(-.5, .5),ylim=c(1, 13))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+   theme(axis.text.x = element_text(size=12))+   theme(axis.title.x = element_text(size=12))+      
  geom_point( colour = 'black', size = 2.75)
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()

variables2<-variables[17:26,]

title<-paste("Effects of Bill Attributes on RPS Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_messaging",".pdf",sep=""), width=6.33, height=3.94) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("4-point Scale for Bill Support")+
  theme(axis.title.y = element_blank())+
  theme(axis.text.y=element_text(size=10, face=rev(variables2$fontface)))  +
   scale_x_continuous(limits = c(-.5, .5),breaks = c(-.4,-.2,0,.2,.4))+
  coord_cartesian(xlim=c(-.5, .5),ylim=c(1, 10))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+   theme(axis.text.x = element_text(size=12))+   theme(axis.title.x = element_text(size=12))+     
  geom_point( colour = 'black', size = 2.75)
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()




#####
### Main Analysis, using binary scale for outcome variable #####
#####

## OLS model without covariates
reg <- lm(billsup_binary~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+ factor(Randomization4b) + factor(Randomization5b), data=master) #+factor(ideo7)
summary(reg)
## OLS model with covariates
reg2 <- glm(billsup_binary~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+ factor(Randomization4b) + factor(Randomization5b), data=master, family = binomial(link = "probit")) #+factor(ideo7)

##  Supplementary Table 5
## comparing models from ols and probit
stargazer(reg,reg2)


## Plot
variables<-read.csv(file="variables.csv")
variables$coef<-NA
variables$SE<-NA

#messenger
variables$coef[24]<-0
variables$coef[25]<-coef(reg)[2]
variables$coef[26]<-coef(reg)[3]

#cost
variables$coef[3]<-0
variables$coef[4]<-coef(reg)[4]
variables$coef[5]<-coef(reg)[5]

#jobs
variables$coef[8]<-0
variables$coef[9]<-coef(reg)[6]
variables$coef[10]<-coef(reg)[7]

#air
variables$coef[13]<-0
variables$coef[14]<-coef(reg)[8]

#climate
variables$coef[18]<-0
variables$coef[19]<-coef(reg)[9]
variables$coef[20]<-coef(reg)[10]
variables$coef[21]<-coef(reg)[11]

SE<-summary(reg)$coefficients[, 2]
#messenger
variables$SE[24]<-0
variables$SE[25]<-SE[2]
variables$SE[26]<-SE[3]

#cost
variables$SE[3]<-0
variables$SE[4]<-SE[4]
variables$SE[5]<-SE[5]

#jobs
variables$SE[8]<-0
variables$SE[9]<-SE[6]
variables$SE[10]<-SE[7]

#air
variables$SE[13]<-0
variables$SE[14]<-SE[8]

#climate
variables$SE[18]<-0
variables$SE[19]<-SE[9]
variables$SE[20]<-SE[10]
variables$SE[21]<-SE[11]
variables<-variables
variables$scale<-"Dichotomous Indicator for Bill Support"

order<-1:26
variables <- transform(variables, Categories=reorder(Categories, -order) ) 

title<-paste("Effects of Bill Attributes on RPS Bill Support", sep="")



variables2<-variables[2:14,]

title<-paste("Effects of Bill Attributes on Republicans' Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_binary_design",".pdf",sep=""), width=4, height=5.25) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("Dichotomous Indicator for Bill Support")+
  # theme(axis.text.y=element_text(size=7, face=rev(variables2$fontface)))  +
   scale_x_continuous(limits = c(-.25, .25),breaks = c(-.2,-.1,0,.1,.2))+
 coord_cartesian(xlim=c(-.25, .25),ylim=c(1, 13))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+  
 theme(axis.title.y = element_blank())+
    theme(axis.text.y = element_blank())+   theme(axis.text.x = element_text(size=12)) +   theme(axis.title.x = element_text(size=12))+    
  geom_point( colour = 'black', size = 2.75)
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()

variables2<-variables[17:26,]

title<-paste("Effects of Bill Attributes on Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_binary_messaging",".pdf",sep=""), width=4, height=3.94) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("Dichotomous Indicator for Bill Support")+
 # theme(axis.text.y=element_text(size=7, face=rev(variables2$fontface)))  + 
 scale_x_continuous(limits = c(-.25, .25),breaks = c(-.2,-.1,0,.1,.2))+
  coord_cartesian(xlim=c(-.25, .25),ylim=c(1, 10))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") +   theme(axis.text.x = element_text(size=12))+   theme(axis.title.x = element_text(size=12))+   
  theme(axis.title.y = element_blank())+
    theme(axis.text.y = element_blank())+
  geom_point( colour = 'black', size = 2.75)
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()
 
 ## both scales
variables2<-rbind(variables1, variables)

fontface<-rev(variables2$fontface[1:52])

order<-1:52

variables2$scale <- factor(variables2$scale,
         levels = c("Bill Support on Four-Point Scale", "Dichotomous Indicator for Bill Support"))


variables2 <- transform(variables2, Categories=reorder(Categories, -order) ) 


title<-paste("Effects of Bill Attributes on Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_bothscales",".pdf",sep="")) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("Difference in Bill Support")+
  theme(axis.title.y = element_blank())+
  theme(axis.text.y=element_text(size=7, face=rev(variables2$fontface)))  +
  coord_cartesian(ylim=c(0, 27))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+  
  geom_point( colour = 'black', size = 2.75)+
   facet_grid( .~scale, scales = "free" )#+
  #        theme(axis.text.y=element_blank())+
  #ggtitle(title)
  
dev.off()




#####
## Supplementary Table 7
## Heterogeneity across regions
#####

##Compare models across RPS and non-RPS states
reg.regions <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b)+factor(educ)+factor(gender)+age+factor(race) +income2+factor(REGION4)  , data=master) # note: includes leaners


##Compare models across RPS and non-RPS states
reg.east <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b)+factor(educ)+factor(gender)+age+factor(race) +income2  , data=master[master$REGION4=="Northeast",]) # note: includes leaners

reg.midwest <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b) +factor(educ)+factor(gender)+age+factor(race) +income2 , data=master[master$REGION4=="Midwest",]) # note: includes leaners

reg.south <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b) +factor(educ)+factor(gender)+age+factor(race) +income2 , data=master[master$REGION4=="South",]) # note: includes leaners

reg.west <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b)+factor(educ)+factor(gender)+age+factor(race) +income2  , data=master[master$REGION4=="West",]) # note: includes leaners

stargazer(reg.east, reg.midwest,reg.south, reg.west)


#####
### Supplementary Table 6
## Heterogeneity by party
#####

## OLS model without covariates for Democrats
reg.dems <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b) , data=master[master$pid3_leaners=="democrat"| master$pid3a=="Democrat",]) # note: includes leaners
summary(reg)

## OLS model with covariates for Democrats
reg.dems.covariates <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b)+factor(educ)+factor(gender)+age+factor(race) +income2 , data=master[master$pid3_leaners=="democrat"| master$pid3a=="Democrat",]) # note: includes leaners

## OLS model without covariates for Republicans
reg.reps<-lm(billsup ~ factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b), data=master[master$pid3_leaners=="republican" | master$pid3a=="Republican",]) # note: includes leaners

## OLS model with covariates for Republicans
reg.reps.covariates <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b)+factor(educ)+factor(gender)+age+factor(race) +income2, data=master[master$pid3_leaners=="republican"| master$pid3a=="Republican",]) # note: includes leaners

stargazer(reg.dems, reg.dems.covariates,reg.reps, reg.reps.covariates)

#######
## Figure 3
########

###### support among D's ####
reg <- lm(billsup~factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b) , data=master[master$pid3_leaners=="democrat"| master$pid3a=="Democrat",]) # note: includes leaners
summary(reg)


variables<-read.csv(file="variables.csv")
variables$coef<-NA
variables$SE<-NA

#messenger
variables$coef[24]<-0
variables$coef[25]<-coef(reg)[2]
variables$coef[26]<-coef(reg)[3]

#cost
variables$coef[3]<-0
variables$coef[4]<-coef(reg)[4]
variables$coef[5]<-coef(reg)[5]

#jobs
variables$coef[8]<-0
variables$coef[9]<-coef(reg)[6]
variables$coef[10]<-coef(reg)[7]

#air
variables$coef[13]<-0
variables$coef[14]<-coef(reg)[8]
  Support among State Legislators   
  #climate
variables$coef[18]<-0
variables$coef[19]<-coef(reg)[9]
variables$coef[20]<-coef(reg)[10]
variables$coef[21]<-coef(reg)[11]

SE<-summary(reg)$coefficients[, 2]
#messenger
variables$SE[24]<-0
variables$SE[25]<-SE[2]
variables$SE[26]<-SE[3]

#cost
variables$SE[3]<-0
variables$SE[4]<-SE[4]
variables$SE[5]<-SE[5]

#jobs
variables$SE[8]<-0
variables$SE[9]<-SE[6]
variables$SE[10]<-SE[7]

#air
variables$SE[13]<-0
variables$SE[14]<-SE[8]

#climate
variables$SE[18]<-0
variables$SE[19]<-SE[9]
variables$SE[20]<-SE[10]
variables$SE[21]<-SE[11]

variables$party<-"Democrats"
variables1<-variables
fontface<-rev(variables$fontface[1:26])

order<-1:26

variables <- transform(variables, Categories=reorder(Categories, -order) ) 
variables2<-variables[2:14,]

title<-paste("Effects of Bill Attributes on Democrats' Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_democrats_design2",".pdf",sep=""), width=6.33, height=5.25) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("4-point Scale for Bill Support")+
  theme(axis.title.y = element_blank())+
  theme(axis.text.y=element_text(size=10, face=rev(variables2$fontface)))  +
  coord_cartesian(xlim=c(-.5, .5),ylim=c(1, 13))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+  
  geom_point( colour = 'black', size = 2.75)+   theme(axis.text.x = element_text(size=12))+   theme(axis.title.x = element_text(size=12))
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()

variables2<-variables[17:26,]

title<-paste("Effects of Bill Attributes on Democrats' Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_democrats_messaging",".pdf",sep=""), width=6.33, height=3.94) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("4-point Scale for Bill Support")+
  theme(axis.title.y = element_blank())+
  theme(axis.text.y=element_text(size=10, face=rev(variables2$fontface)))  +
  coord_cartesian(xlim=c(-.5, .5),ylim=c(1, 10))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+  
  geom_point( colour = 'black', size = 2.75)+   theme(axis.text.x = element_text(size=12))+   theme(axis.title.x = element_text(size=12))
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()

#### support among R's #####
reg<-lm(billsup ~ factor(Randomization1b)+factor(Randomization2b)+factor(Randomization3b)+factor(Randomization4b)+factor(Randomization5b), data=master[master$pid3_leaners=="republican" | master$pid3a=="Republican",]) # note: includes leaners
summary(reg)

variables<-read.csv(file="variables.csv")
variables$coef<-NA
variables$SE<-NA

#messenger
variables$coef[24]<-0
variables$coef[25]<-coef(reg)[2]
variables$coef[26]<-coef(reg)[3]

#cost
variables$coef[3]<-0
variables$coef[4]<-coef(reg)[4]
variables$coef[5]<-coef(reg)[5]

#jobs
variables$coef[8]<-0
variables$coef[9]<-coef(reg)[6]
variables$coef[10]<-coef(reg)[7]

#air
variables$coef[13]<-0
variables$coef[14]<-coef(reg)[8]

#climate
variables$coef[18]<-0
variables$coef[19]<-coef(reg)[9]
variables$coef[20]<-coef(reg)[10]
variables$coef[21]<-coef(reg)[11]

SE<-summary(reg)$coefficients[, 2]
#messenger
variables$SE[24]<-0
variables$SE[25]<-SE[2]
variables$SE[26]<-SE[3]

#cost
variables$SE[3]<-0
variables$SE[4]<-SE[4]
variables$SE[5]<-SE[5]

#jobs
variables$SE[8]<-0
variables$SE[9]<-SE[6]
variables$SE[10]<-SE[7]

#air
variables$SE[13]<-0
variables$SE[14]<-SE[8]

#climate
variables$SE[18]<-0
variables$SE[19]<-SE[9]
variables$SE[20]<-SE[10]
variables$SE[21]<-SE[11]

fontface<-rev(variables$fontface[1:26])

order<-1:26
variables <- transform(variables, Categories=reorder(Categories, -order) ) 

variables2<-variables[2:14,]

title<-paste("Effects of Bill Attributes on Republicans' Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_republicans_design",".pdf",sep=""), width=4, height=5.25) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("4-point Scale for Bill Support")+
  # theme(axis.text.y=element_text(size=7, face=rev(variables2$fontface)))  + scale_x_continuous(limits = c(-.5, .5),breaks = c(-.5,-.25,0,.25,.5))+
 coord_cartesian(xlim=c(-.5, .5),ylim=c(1, 13))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+  
 theme(axis.title.y = element_blank())+
    theme(axis.text.y = element_blank())+
  geom_point( colour = 'black', size = 2.75)+   theme(axis.text.x = element_text(size=12))+   theme(axis.title.x = element_text(size=12))
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()

variables2<-variables[17:26,]

title<-paste("Effects of Bill Attributes on Republicans' Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_republicans_messaging",".pdf",sep=""), width=4, height=3.94) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("4-point Scale for Bill Support")+
 # theme(axis.text.y=element_text(size=7, face=rev(variables2$fontface)))  + scale_x_continuous(limits = c(-.5, .5),breaks = c(-.5,-.25,0,.25,.5))+
  coord_cartesian(xlim=c(-.5, .5),ylim=c(1, 10))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + 
  theme(axis.title.y = element_blank())+
    theme(axis.text.y = element_blank())+
  geom_point( colour = 'black', size = 2.75)+   theme(axis.text.x = element_text(size=12))+   theme(axis.title.x = element_text(size=12))
  #        theme(axis.text.y=element_blank())+
 #+ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
 dev.off()


### merged figure for both parties
variables$party<-"Republicans"

variables2<-rbind(variables1, variables)

fontface<-rev(variables2$fontface[1:52])

order<-1:52



variables2 <- transform(variables2, Categories=reorder(Categories, -order) ) 


title<-paste("Effects of Bill Attributes on Bill Support", sep="")
library(ggplot2)
pdf(paste("RPS_Bill_regression_results_bothparties",".pdf",sep="")) # for paper
ggplot(data = variables2, aes(x = coef, y = Categories)) +
  geom_errorbarh(data = variables2, aes(y = Categories, 
                                       xmin = coef - 1.96*SE, xmax = coef + 1.96*SE, height=.25),
                 colour = 'black')+ xlab("4-point Scale for Bill Support")+
  theme(axis.title.y = element_blank())+
  theme(axis.text.y=element_text(size=7, face=rev(variables2$fontface)))  +
  coord_cartesian(xlim=c(-.5, .5),ylim=c(0, 27))+ geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") + theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black"))+  
  geom_point( colour = 'black', size = 2.75)+
   facet_grid( .~party )+
  #        theme(axis.text.y=element_blank())+
 ggtitle(title)+theme(plot.title = element_text(hjust = 0.5))
dev.off()

